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ABSTRACT 

We study a slightly rotating accretion flow onto a black hole, using the fully three dimen- 
sional (3-D) numerical simulations. We consider hydrodynamics of an inviscid flow, assuming a 
spherically symmetric density distribution at the outer boundary and a smafl, latitude-dependent 
angular momentum. We investigate the role of the adiabatic index and gas temperature, and 
the flow behaviour due to non-axisymmetric effects. Our 3-D simulations conflrm axisymmetric 
results: the material that has too much angular momentum to be accreted forms a thick torus 
near the equator and the mass accretion rate is lower than the Bondi rate. In our previous 
study of the 3-D accretion flows, for 7 = 5/3, we found that the inner torus precessed, even for 
axisymmetric conditions at large radii. The present study shows that the inner torus precesses 
also for other values of the adiabatic index: 7 = 4/3, 1.2, and 1.01. However, the time for the 
precession to set increases with decreasing 7. In particular, for 7 = 1.01 we flnd that depending 
on the outer boundary conditions, the torus may shrink substantially due to the strong inflow of 
the non- rotating matter and the precession will have insuflicient time to develop. On the other 
hand, if the torus is supplied by the continuous inflow of the rotating material from the outer 
radii, its inner parts will eventually tilt and precess, as it was for the larger 7's. 

Subject headings: accretion, accretion discs - black hole physics - galaxies: active 



1. Introduction 

The black holes capturing the gas from their 
vicinity are ubiquitous in the Universe, and accre- 
tion onto the central black hole is the source of 
energy that powers the brightest sources on our 
sky, from quasars, through the gamma ray bursts, 
to the X-ray binaries. The fast rotation of ma- 
terial before it becomes captured by the central 
object leads to the formation of a thin accretion 
disk, where the gravitational energy is effectively 
dissipated through the viscous stresses and radi- 
ated away in the form of high energy radiation. 
On the other hand, if the gas does not rotate. 
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or rotates very slowly, the viscous disk does not 
form and the material falls radially onto the cen- 
ter. The rate of accretion is then equal to the 
Bondi accretion rate, derived analytically (Bondi 
1952). The flows with a small, latitude dependent 
angular momentum were studied numerically e.g. 
in Proga & Begelman (2003a,b) in 2D, and Janiuk 
et al. (2008) in 3D, and showed that in this case 
the accretion rate is on the order of 20-30% of the 
Bondi rate. The astrophysical situation that may 
be described by such solutions is the radiatively 
inefficient ffow of gas captured by the supermas- 
sive black holes in low luminosity AGNs, or the 
quiescent state of an X-ray binary. 

The simulations show, that if the torus is 
formed due to rotation, the gas cannot accrete 
unless some angular momentum transport mech- 
anism is introduced. Therefore no accretion pro- 
ceeds close to and through the equatorial plane, 
and the material accumulated in this region either 
flows out due to the centrifugal forces and gas 
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pressure, or eventually turns to smaller latitudes 
and flows radially through the polar region. The 
recent three dimensional studies showed, that in 
the latter case the amount of gas turning in the 
direction above or below the torus may depend on 
the azimuth, and in consequence lead to the tilt 
of this torus and its precession. 

The present work is an extension of the previ- 
ous study (Janiuk et al. 2008) to check how the 
gas properties, described by its equation of state, 
will influence the results. The general form of the 
EOS assumed here is the polytropic relation, and 
the index of 7 = 5/3, appropriate for gas pressure 
dominated system, is well suited for example in 
the inactive AGNs. On the other hand, the GRBs, 
which result from accretion of the massive stellar 
envelope, being a mixture of nucleons, neutrinos, 
electron-positron pairs and trapped photons (e.g. 
Janiuk et al. 2004; 2007), will rather be better de- 
scribed by a relativistic EOS with 7 = 4/3, than 
the ideal gas EOS. The still smaller values of 7, 
discussed in Moscibrodzka & Proga (2008), can 
account for the nearly isothermal gas, suggested 
to constitute the protogalactic disks. 

Therefore we consider here the 3-dimensional 
hydrodynamical models with a range of 7 in- 
dices describing different gas microphysics. We 
concentrate on the initially axisymmetric condi- 
tions, which were shown to produce the non- 
axisymmetric results in the process of time evo- 
lution. The content of the article is as follows. In 
Section [2j we describe the method used in our cal- 
culations, to determine the initial conditions and 
subsequent evolution of the flow. In Section [3j we 
present the simulations results for the 3-D models 
with various adiabatic indices (Sec. 3.1). We dis- 
cuss the torus precession and the influence of the 
outer boundary conditions (Sec. |3.2[ ), which oc- 
cured to be important, especially for models with 
smaller 7. The models with various gas temper- 
atures are presented in Sec. 3.3 We discuss our 
results in Section [H 

2. Method 

2.1. Initial conditions 

In the analytical formula of Bondi (1952), the 
steady state solution to the equations of mass and 
angular momentum conservation is parameterized 
with the two quantities at inflnity: density, p^^^ 



and sound speed, c^^oo- The solution is given by 
the accretion rate of: 



Mb = A47ri?|pooCs 



(1) 



where the dimensionless parameter A depends only 
on the adiabatic index: 
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for 1 < 7 < 5/3. In the limit of 7 ^ 1.0, this 
parameter approaches: 



A= Jexp(3/2), 



(3) 



while in the limit of 7 = 5/3 the accretion rate 
is calculated from the transonic solution in the 
pseudonewtonian potential, PW, (Paczyhski & 
Wiita 1980), and depends on the sound speed 
at inflnity (see Proga & Begelman 2003a for the 
derivation from the density and sound speed at 
the sonic radius): 



The Bondi radius is equal to 
GM 

— -7^ . 



(4) 
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The transonic solution is obtained for 7 < 5/3 
in the Newtonian potential, and the sonic radius 
is located at: 

GM 

with the sound speed in the sonic radius equal to: 

(7) 



Cs('^s) — Cg^c 
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The above relation gives no sonic point for 7 = 
5/3, and in this case the transonic solution is found 
for the PW potential, with the sonic radius equal 
to: 

- ^ + [2 GM 

We model the initial conditions for the spher- 
ically symmetric Bondi inflow, solving iteratively 
the Bernoulli equation in 1-D: 
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where 



B = H 
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GM 



(9) 



(10) 



is the PW potential, Rs = 2GM/c^ is the 
Schwarzschild radius, and 
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is the enthalpy. 

We parameterize our model with poo and Cs,oo- 
The latter is conveniently used to define a dimen- 



sionless radius, 



Rs/Rb 



2clJc\ De- 



pending on 7, we calculate the dimensionless ac- 
cretion rate from Eq. |2) Eq. [3) or Eq. [4| Then, 
we solve for the density profile, while the velocity 
is found from the continuity equation: 



Mb 



-47rr v^p. 



(12) 



In the initial condition, we also impose a slight 
rotation of the outermost parts of the Bondi flow. 
This is done by means of the latitude-dependent 
angular momentum. 



/ = /o(l - |cos(9|) 



(13) 



where /q is a dimensionless parameter, which cor- 
responds to the circularization radius Rq (i.e. the 
radius at which GM/r^ = v'^/r). Therefore, 
our initial velocity vq is everywhere equal to 0.0, 
whereas the velocity 'U^ is initially non-zero only in 
a fixed, quasi-conical region, with its radial exten- 
sion limited by the sound speed at infinity. This 
region is defined as follows: 
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for 
for 



I'^r I ^ Cg^QQ 
\Vr\ < Cs.oo 

(14) 



2.2. Boundary conditions 

The boundary conditions are set at the outer 
radius of the computational domain. We impose 
a spherically symmetric outer boundary condition 
for the fiow density, ^(rout), and energy density, 
^(rout), to be equal to these quantities determined 
from the initial analytical soulution for the Bondi 
fiow. Because the outer radius of our grid is fi- 
nite, the density at the outer edge is not exactly 



equal to poo, but roughly twice larger. Neverthe- 
less, the outer radius of 1.2 Rb is suflicient and 
such a condition matches well the outer, infinite 
and stationary sphere of gas, for which all the vari- 
ables are determined analytically, with the inner, 
numerically simulated region. 

The boundary conditions for the velocity field 
are not specified otherwise than these incorpo- 
rated in the ZEUS-MP code. We use the free in- 
fiow/outfiow condition for the radial velocity, the 
refiection symmetry with respect to the polar axis, 
and periodic boundary condition in the azimuthal 
direction. 

Note, that in this way in every time step we up- 
date the value of density at the outer radius, but 
we do not update the specific angular momentum. 
In other words, the physical situation that we con- 
sider, is a spherically symmetric, stationary cloud 
of gas accreting from the infinity (where the den- 
sity and temperature are given by poo and Cg^oo) 
onto the central massive black hole. This cloud 
was at time t = perturbed by imposing some 
small amount of angular momentum, which made 
the cloud evolve. 

2.3. Time evolution of the flow 

In our calculations we use the 3-D code ZEUS- 
MP (Stone & Norman 1992; Hayes & Norman 
2003), modified by ourselves to incorporate the 
PW gravitational potential. The computations 
were performed on the multi-processor computer 
cluster machines (see Catlett et al. 2007). The 
ZEUS-MP code solves the equations of hydrody- 
namics: 



dp 
Itt 
dv 

p4(-) + i^Vv = 
at p 



■ pVv = 



-VP + 



(15) 



(16) 
(17) 



where p is the gas density, e is the internal energy 
density, P = (7 — l)e is the gas pressure, and v is 
the velocity of the fiow. 

The other model parameters, Cg^oo, Poo, ^bh 
are chosen such that the ratio of the Bondi radius 
to the Schwarzschild radius is fixed, and is equal 
to 1000, 300 or 100, depending on the model. The 
radial grid is in the range from rin = l.bRs^ to 
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^out = 1.2i?B- The rotation of the flow is param- 
eterized by the value of Rq, and for most of the 
models it is equal to 0.1 (in the units of Bondi 
radius). 

We use the spherical coordinate system, RTP, 
and the resolution in r-direction was 140 zones, 
with dr\j^i/dr\ = 1.05, in ^-direction it was 96 
zones, and in (/)-direction it was 32 zones, with 
dOj-^i/dOj = d(j)\^j^i / d(j)\^ = 1.0. The range of the 
grid in and (j) directions is from to tt and from 
to 27r, respectively. 

3. Results 

The hydrodynamical computations of an ax- 
isymmetric, 2.5-D model of an accretion flow with 
low angular momentum were presented in PB03, 
and continued in Moscibrodzka & Proga (2008). 
In the 3-D analysis presented in Janiuk, Proga 
& Kurosawa (2008), we recalculated their 2.5-D 
models for one chosen value of the adiabatic in- 
dex, 7 = 5/3, and the sound speed at inflnity cor- 
responding to the Bondi radius equal to 1000 Rs. 
In that paper, we confirmed that the 3-D effects, 
such as the nonaxisymmetric distribution of the 
angular momentum in the accreting fluid, play an 
important role for a resulting rate of mass accre- 
tion with respect to the Bondi one. We also found 
that a non-axisymmetric torus is subject to the 
tilt and precession due to the acoustic instabilities 
in the innermost gas. These instabilities appear 
when any small asymmetry in the angular mo- 
mentum distribution arises during the torus evo- 
lution. In the dynamical timescale, they lead to 
the torus misplacement from the equatorial plane 
and its subsequent precession. We showed that 
the instability develops in a highly supersonic flow, 
whereas for small Mach numbers it is suppressed. 
We attributed this behaviour with the Papaloizou 
& Pringle (1995) type of instability, which for low 
order modes is driven by the Kelvin-Helmholtz 
mechanism. 

In all of our models, the simulations start from 
a spherically symmetric gas cloud around a black 
hole, with density and velocity distributions de- 
rived from the Bondi solution. The matter lo- 
cated far from the black hole contains speciflc an- 
gular momentum that exceeds the critical value, 
/crit = 2i?sc, at the equatorial plane and is de- 



creasing towards the polar regions (cf. Eq. 14). 



The time evolution of the system with such ini- 
tial conditions proceeds flrst through a short phase 
of the purely radial Bondi accretion, which du- 
ration depends on the assumed distance of the 
initially rotating gas from the central black hole. 
Then, the evolution of the flow switches to a long- 
term phase of torus accretion. In this phase, the 
gas settled near the equatorial plane is supported 
against gravity by the gas pressure and rotation, 
and the rate of accretion on the black hole, Min, 
decreases below ~ 30% of the Bondi rate. The 
reason is that the material from the polar regions 
is still accreting radially, while the torus mate- 
rial is mainly rotating, and is either outflowing, 
or accretes after turning slighlty towards one of 
the poles. 

Most of our simulation runs lasted up to about 
1.5 X 10^ dynamical time at the inner radius, t = 
tdyn(nn) = '^'7^ / {^in) ^ whcrc Qk is the Keplcriau 
velocity. 

The axisymmetry of the solution breaks after a 
certain time, depending on the model parameters 
(see below). When this happens, the inner torus 
becomes tilted with respect to the equator and it 
starts precessing. This may also be accompanied 
by the fluctuations of the accretion rate. 

3.1. Evolution of the flow for various adi- 
abatic indices 

We performed the runs for the axisymmetric 
initial conditions for a range of 7 and gas temper- 
ature (coo, reflected by the value of Rb/Rs)- The 
models are summarized in Table [3J1 

In Figure [l] we show the time evolution of the 
accretion rate through the inner boundary, Min (in 
units of the Bondi accretion rate), for various adi- 
abatic indices. As the Figure shows, the accretion 
rate onto BH varies in time. Before the rotating 
material approaches the black hole, Min is equal 
to the Bondi accretion rate for all models. Once 
the gas starts rotating also in the innermost parts, 
the accretion rate drops to a small fraction of the 
Bondi rate. The moment when this happens, i.e. 
the end of the transient phase of the purely radial 
accretion, depends on the adiabatic index, and is 
in the range from about t = 3.2 x 10^ for 7 = 5/3, 
up to about t = 9.5 x 10^ for 7 = 1.01. 

In all the models, the accretion rate after the 
torus has formed, is on average low, and does not 
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Mo 


7 


Rb/Rs 


^end 


M 


Prec 


A 


1.01 


1000 


1.1 X 10"^ 


0.12-1.55 




Av 


1.01 


1000 


1.2 X 10"^ 


0.1-0.75 


yes 


B 


1.2 


1000 


1.25 X 10"^ 


0.09-0.61 


yes 


C 


4/3 


1000 


1.25 X 10"^ 


0.10-0.28 


yes 


D 


5/3 


1000 


1.25 X 10"^ 


0.18-0.25 


yes 


E 


5/3 


300 


8.1 X 10^ 


0.30-0.42 


no 


F 


5/3 


100 


6.2 X 10^ 


0.34-0.62 


no 


G 


4/3 


300 


6.2 X 10^ 


0.12-0.32 


no 


H 


4/3 


100 


6.1 X 10^ 


0.05-0.30 


no 



Table 1: Summary of the evolutionary models. 
The time is given in units of dynamical time at 
Tin- The accretion rate M is given in the units of 
Mb. 
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Fig. 1. — Time evolution of the mass accretion 
rate through the inner boundary (Min), in units 
of the Bondi rate (Mb). The models are: A (red 
line), B (green line), C (black line), D (cyan line). 



exceed 30% of the Bondi rate. The fluctuations of 
the accretion rate for 7 = 1.2, 4/3 and 5/3 are not 
signiflcant: they are on the order of 0.02 — 0.04Mb. 
However at the end of the simulation, for 7 = 1.2 
the accretion rate approaches and exceeds 0.5Mb. 
For the almost isothermal model (7 = 1.01), this 
happens much earlier. The small accretion rate is 
kept only for a very short time, t ^ 3.9 x 10^, and 
after that Min starts rising, to reach and exceed 
the Bondi rate. This is accompanied by very large 
fluctuations of the accretion rate, in a form of pro- 
nounced flares with amplitude of 0.3 — 0.5Mb. 

In the Figure |2j we plot the evolution of the 
angular momentum flux through the inner bound- 
ary: Lin = / IspecP'^rds^ in uuits of the critical 
angular momentum /crit and renormalized by the 
value of the Bondi accretion rate. At the begin- 
ning of the simulation, Lin = 0, because in the 
vicinity of the black hole the gas is not rotating. 
Once the rotating matter reaches the inner bound- 
ary, the angular momentum starts accreting to the 
center, and Lin sharply increases. However, after 
several orbital cycles the outflow begins and the 
net radial velocity as well as the density near the 
polar regions drop, so Lin decreases to a moderate 
value of ~ 0.1 (in the units of /cot Mb). This is 
the case for the gas and ratiation pressure domi- 
nated models, 7 = 5/3 and 7 = 4/3, and also for 
the model with 7 = 1.2 for most of the simulation. 
Only for the model with 7 = 1.01, the angular mo- 
mentum flux rises after t ^ 3.9 x 10^ and strongly 
oscillates. 

Note that the quantity plotted in the Figure [2] is 
a flux of total angular momentum, not specific an- 
gular momentum. This corresponds to the amount 
of angular momentum which may be transferred 
to the black hole and used to spin it up. Fig- 
ure [2] can be used to calculate the total angu- 
lar momentum which the black hole could gain 
during our simulation. We find, after the proper 
unit conversions, that this number is extremely 
low: a = {cJ)/GM'^ ^ 4 x 10"^ - 10-^ where 
J = Lin At. For the isothermal model, 7 = 1.01, 
the amount of the total angular momentum ac- 
creted onto the black hole will be 2-3 times larger 
(cf. Fig. [1]). Still, it does not lead to the spin- up of 
the black hole, as the material which is being ac- 
creted is the very slowly rotating gas mostly from 
the polar regions. In the models with 7 = 5/3, 4/3 
and 1.2, the accretion rate and Lin are only slightly 
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variable. 

To understand the behaviour of the accretion 
rate Min, we need to study the structure of the 
flow. In Figures [3| [5| [7| and |9| we plot the contour 
maps of the density and velocity fields in the in- 
nermost regions of the fiow (i.e. up to 20 Rs), for 
the models A,B,C and D, respectively. The maps 
are plotted in three time-snapshots, and in the two 
views: equatorial plane and side-on. 

In Figure [3] we show the density contour maps 
for the inner region of the accretion fiow, for the 
model A (7 = 1.01). The maps also show the 
velocity field, and are plotted for several snapshots 
during the evolution: t ^ 3 x 10^, t ^ 6.3 x 10^ 
and t ^ l.lx 10^. The panels on the left, show the 
side-on view of the inner torus, i.e. the x — z-plane, 
while the panels on the right show the top-view, 
i.e. the equatorial plane. 

The torus, which formed after the transient 
phase of the Bondi accretion, is very dense: the 
maximum density in the equatorial plane exceeds 
8 X lO^poo. At the time t - 3 x 10^, the po- 
lar regions are clean, and most of the material 
is settled near the equator. The radial infall of 
material onto the central black hole goes mainly 
through the poles, while the torus material is tur- 
bulent and due to the centrifugal barrier the gas 
may fiow outwards. At this time, the rate of accre- 
tion through the inner boundary, Min, is the low- 
est and is below 20% of the Bondi rate. The rapid 
variations of the accretion rate are caused by the 
turbulences. Later, the variations are enhanced 
by the azimuthal perturbations, which grow in the 
density field in the innermost region of the torus 
(see the middle-right panel of the Figure [sj for 
time t - 6.3 X 10^). As was shown in Figure 
[1] in the late phase of the evolution in model A, 
the accretion rate Min steeply rises, reaches the 
Bondi rate, and rapidly varies about this value. 
However, the fiow does not return completely to a 
spherically symmetric configuration of the Bondi 
type. The inner torus still exists, but becomes 
much smaller and its density is lower. The maxi- 
mum density in the equatorial plane is now about 
^ 2 X lO^Poo, while the density in the polar regions 
becomes larger. 

Now the polar material forms a kind of streams, 
in which the density is about 100 times larger in 
the end of the run than it was during the phase 
of 'clean poles'. This large density, in addition to 



the highly supersonic radial velocities (see Fig. [4] 
below), are the cause for the net accretion rate 
to be large again. The streams are formed by 
the material that has too large angular momen- 
tum to accrete radially, but fiows laminarly over 
the torus and turns toward the poles. The 'arcs' 
visible on the density map are the results of the 
obligue shocks, where the vq velocity component 
is changing. 

We note, that such kind of behaviour, i.e. the 
drammatic rise of the accretion rate in the nearly 
isothermal mode, 7 = 1.01, was not found in the 
study by Moscibrodzka & Proga (2008), however, 
their 2-D runs were performed within twice longer 
timescale, with Tend ~ 4 x 10^. In their simula- 
tions, a twice smaller rotation parameter was used, 
lo = 0.05, so that the accretion rate Min in their 
models did not drop below ~ 0.4Mb. We tested 
that the difference in the rotation parameter is not 
crucial: we made the test calculations with other 
values of Iq and qualitatively the results were sim- 
ilar. Also, we checked that it is not the role of the 
3-D effects, that make the nearly isothermal torus 
smaller and weaker after a short time in our simu- 
lation. The same model was tested in the 2.5-D ax- 
isymmetric configuration with ZEUS-MP, and we 
found that the mean accretion rate also increased 
to about Mb after t ~ 6 x 10^, with the fiaring be- 
haviour being even more pronounced than in our 
3-D case. 

What we found crucial for the behaviour of the 
fiow at small 7, was the treatment of the outer 
boundary condition for the azimuthal velocity. In 
their simulations with ZEUS 2-D, Moscibrodzka & 
Proga (2008) imposed the update of specific angu- 
lar momentum at the outer boundary in every time 
step. In our simulations, we imposed the rotation 
of the outer regions only in the initial condition, 
to mimic the transient cloud of gas, that contains 
some angular momentum, passing by the vicinity 
of the galaxy center. When the cloud sinks into 
the Bondi sphere, the remaining gas at infinity 
no longer rotates. The Bondi sphere is still mod- 
eled with the boundary condition for density, and 
matched with the analytical Bondi solution. 

We note here, that despite the rise of the ac- 
cretion rate to the value of the Bondi rate, the 
torus in the innermost parts of the fiow does not 
completely disappear. The torus shrinks in size 
and is less dense, because of the massive infiow of 
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the non-rotating material from the outer bound- 
ary. However, it stih acts as an obstacle for the 
gas flowing radially onto the black hole, which re- 
sults in the rapid fluctuations of the accretion rate 
at the inner boundary. 

In Figure |4] we show the contour plots of the 
total Mach number, Mtot = |'^|/cs (left) and the 
radial Mach number, = v^/c^ (right), for 
7 = 1.01. The color contours show the constant 
Mach number of 0.5, 1.0 and 1.5 (red, green and 
blue, respectively), while the grey colors show the 
three chosen constant denity contours, to mark the 
position of the inner torus. In model A, which is 
nearly isothermal, most of the flow is highly su- 
personic. Starting from the Bondi solution, the 
sonic radius is located at a very large distance, 
rs = 0.49i?B, and the total Mach number is every- 
where by deflnition equal to the radial one. When 
the torus forms, it is supported by rotation, and 
the radial velocity in the torus is smaller than the 
sound speed. Only in the innermost radii, below 

2.6 X 10"^i?B = 2.6i?s, the radial Mach number 
is greater than 1.0. The sonic radius for the total 
Mach number is however at the same place as be- 
fore, and the net velocity of the flow is supersonic. 
At the time t ~ 6.3 x 10^ (middle panels of the 
Figure [4| , the sonic radius in the equatorial re- 
gion shrinks. The radial and the total Mach num- 
bers are both greater than 1.0 close to the poles. 
The supersonic velocity components closer to the 
equatorial plane are mostly non-radial. Later dur- 
ing the simulation, this trend reverses. The sonic 
surface for Mtot closes again at ~ 0.5i?B, and 
it coincides with the sonic surface for (not 
shown in the bottom-right flgure, because of the 
zooming-in). Close to the equatorial plane, the re- 
gion with subsonic radial velocities shrinks several 
times, and at the outer edge of the torus, the radial 
velocity is equal to the sound speed. What is also 
worth noticing in Fig. [4j is that the actual size of 
the constant density contours, that has shrunken 
about 10 times in the last time snapshot. 

In Figure |5] we show the flow topology for the 
model B (7 = 1.2). In this model, the density 
of the torus is much lower. Initially, the flow is 
axially symmetric, but at the time t' ^ 6.3 x 10^, 
the axial symmetry is broken. This is because the 
material at one side of the torus (i.e. the azimuth 
(j) ~ 0°), tries to reach the black hole and flows 
in above the equatorial plane, while at the other 



side (i.e. the azimuth (j) ~ 180°), the gas tries to 
flow in below the equatorial plane. This is shown 
by the directions of the arrows that represent the 
velocity fleld. As was already shown in Janiuk et 
al. (2008), such a behaviour leads to the tilt of 
the torus, and here the tilt is clearly visible at the 
time t' - 1.1 X 10^. 

In Figure [6] we show the Mach number contours 
for the model 5, 7 = 1.2. The color contours 
show the constant Mach number of 0.5, 1.0 and 
1.5 (red, green and blue, respectively), while the 
grey colors show the constant denity contours, to 
mark the position of the inner torus. The shape of 
the sonic surface is very different from the previ- 
ous case (in model A), and resembles in shape the 
number "eight". The radial velocity exceeds the 
speed of sound in the polar regions, while near the 
equator, somewhat irregular contour of the total 
Mach number shows the contribution from the su- 
personic, non-radial velocity components. When 
the torus gets tilted, the innermost part of the 
sonic surface tilts as well, and the size of this sur- 
face slightly increases. 

Similar results are obtained for the adiabatic in- 
dices 7 = 4/3 and 7 = 5/3. The density maps are 
shown in Figures [7| and |9j while the Mach number 
contours are shown in Figures [S] and [To) We note 
that qualitatively, the basic pattern of the torus 
evolution is uniform for 7 = 1.2,4/3 and 5/3 in 
the time scale of our simulations. We note that 
there is a systematic decrease in the torus densi- 
ties with increasing 7 (in Fig. [9] the color scale 
had to be changed to make the density maps more 
contrasted). 

We also note here that in the late time evo- 
lution of model seen in Fig. [9] (bottom- left 
panel) , the flow pattern is strongly affected by the 
poles. This is the result of a discontinuity at the 
z-axis, characteristic for the RTF coordinate sys- 
tem. The boundary conditions in the theta direc- 
tion, which we use in our calculations (reflection 
with inversion in 2nd and 3rd velocity component), 
are appropriate as long as the flow is basically ax- 
isymmetric. Because in the late time evolution 
the torus tilts, the z-axis is no longer the symme- 
try axis. Especially, for the model D the tilt is 



largest (see below. Sec. 3.2), so the pole effects for 



this 7 = 5/3 are the strongest at the end of the 
evolution. 

To summarize and compare the evolution of 
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the flow for various adiabatic indices, we show 
the large scale poloidal velocity flelds, plotted for 
models A,B,C and D, at three representative time 
snapshots. Figure 11 shows, from top to bottom, 
the velocity flelds at time t = 3x 10^, 6.3 x 10^ and 
1.1 X 10^, taken at the azimuth ^ = 0. Because 
the model is not axially symmetric, the snapshots 
taken at different azimuths will not look exactly 
the same. Nevertheless, the global pattern, which 
depends mainly on 7, is very similar. The flrst 
column shows the model with 7 = 1.01. Here, 
the torus which formed initially in the equatorial 
region, is accompanied by the departure of the ve- 
locity fleld from the purely radial. The gas turns 
back due to the centrifugal force, however no net 
outflow occurs, and the largest size of the turbu- 
lent region is about 200 i^schw In the last time 
snapshot, for 7 = 1.01 the torus is indeed very 
small, and the size of the turbulent region is only 
about 20 i^schw, hardly visible in the scale of the 
Figure. 

For other models, the turbulent region is much 
larger, and grows with 7. The shock front be- 
tween the outflowing and inflowing gas is clearly 
visible, and this front propagates outwards with 
time. The large outflow of material occurs from 
the equatorial region, and is accompanied by the 
large scale turbulences. 

3.2. Torus precession 

The precession of the torus is discussed below 
in terms of the tilt and twist angles. These angles 
are deflned as (see e.g. Fragile et al. 2007; Janiuk 
et al. 2008): 



/3{r^t) = arccos (^"^^ 



and 



j{r^t) = arccos ( — ^ 



(18) 



(19) 



Initially, the flow is axisymmetric and the only 
component of the angular momentum vector is 
(which sign depends on the direction of the flow ro- 
tation), while Lx and Ly are basically zero. There- 
fore the tilt angle is vanishingly small and the disk 
does not precess. After some time, the non-zero 
Lx and Ly components appear in the flow, while 
|Lz| decreases, and the rotation axis tilts towards 



the X — y plane. Then, Ly and L^ are chang- 
ing periodically, with L^ being almost constant, 
i.e. the rotation axis moves clockwise or counter- 
clockwise, depending on the model. 

The evolution of the tilt and twist angles pro- 
ceeds as follows. The tilt (i.e. the angle between 
the angular momentum vector of the gas and the 
z axis), is initially equal to zero for all radii. Later 
during the simulation the tilt rises strongly in the 
inner parts of the flow, while it is negligible in the 
outer parts. The twist angle is deflned as a cu- 
mulative angle by which the angular momentum 
vector revolves in the x — y plane by the time t. 
Before the disk was tilted, it did not precess, and 
by deflnition the twist was zero everywhere. When 
the tilt increases, the twist angle rises fast in the 
innermost parts of the flow, which we show in the 
Figure 12 (as well as in the following Figures 13, 
14 and 15) as the twist averaged over the radius. 



from Tin to 0.05i?B- Initially, the strongest rise for 
the innermost radii implies the differential preces- 
sion of the torus. The outermost parts of the flow 
have a negligible precession, and above ~ 0.5i^B 
the twist oscillates around zero. The total twist, 
i.e. the average over the radius from rin to rout, is 
smaller than that for the innermost torus. 



In the Figure [12] we show the tilt and twist evo- 
lution for the model with 7 = 5/3. The tilt is 
shown as a function of radius, for several time- 
snapshots. The initial tilt of 180° means that the 
angular momentum vector is not tilted and the 
flow rotates with a negative azimuthal velocity. 
The twist angle is shown as a function of time, 
and the averaging over radius was calculated in 
two ranges: from r = rin to r = 0.05i?B (solid 
line) and to r = 1.2R^ (dashed line). The twist 
angle decreasing from ~ 60° to ~ —80° (inner- 
most torus) as well as from ~ 40° to ~ —40° (to- 
tal) reflects the clockwise precession. Initially, the 
precession is differential, i.e. only the innermost 
parts of the torus tilt and twist. Later, during the 
time evolution, also the outer parts of the flow tilt, 
and the last two time snapshots in the left panel 



of the Figure 12 show that the tilt increases uni- 
formly for all radii. In the right panel of the same 
Figure, the slopes of the dashed and solid lines at 
late times are nearly the same, meaning that after 
the period of a differential precession, the torus is 
precessing as a nearly solid body. 

In Figure [13] we plot the tilt and twist angles for 
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the model with 7 = 4/3. Here again, the initial 
tilt is of 180°, meaning no tilt and negative az- 
imuthal velocity. The twist angle increasing from 
^ —120° to ~ —40° (innermost torus) as well as 
from ~ —70° to ~ 0° (total) shows the counter- 
clockwise, initially differential, and at late times a 
nearly solid body precession. 

In Figure [14] we plot the tilt and twist evolution 
for 7 = 1.2. The precession starts even later than 
in case of 7 = 4/3, and proceeds also counter- 
clockwise. 

For the model with 7 = 1.01, we cannot confirm 
neither tilt no precession of the torus. In Figure 



15, there is no trend of increase for the tilt and 
twist angles, and the large scatter in the two angles 
is caused by the virtually very small components 
of the angular momentum vector, and Ly. We 
find that the tilt fluctuates periodically and no 
clear signs of precession towards any direction can 
be detected. 

We estimated the times when precession started, 
and periods of precession, for models B,C and D. 
The start times are: t = 4 x 10^, 8.5 x 10^ and 
9.5 X 10^, for 7 = 5/3,4/3 and 1.2, respectively. 
The delays in the torus tilt are corresponding to 
the similar delays in the torus formation, in case 
of smaller 7 with respect to 7 = 5/3. On the other 
hand, the precession period can be estimated at 
2.4 X 10"^, 1.6 X 10^ and 6.4 x 10^. Therefore we 
conclude, that the larger the adiabatic index, the 
earlier the torus forms and starts precessing, while 
the longer the precession period. 

For 7 = 1.01 no such estimates could be made, 
because the torus shrank too much before it could 
start precessing. However, we note that on the 
density map (Fig. [sj middle right panel) there is 
a non-axisymmetric structure, which develops in 
the innermost region. This might be a hint for 
an azimuthal (m = 2) instablity mode that starts 
developing at this time in the simulations. There- 
fore we made an additional test calculation, model 
Av^ in which we changed the outer boundary con- 
dition for the velocity and we continuously added 
the angular momentum at the outer radius. In this 
way, the torus could be supported rotationally for 
a very long time during the simulation, and did 
not shrink like in model A. 

The accretion rate dependence on time for this 
model is shown in Figure 16 There is clearly no 



rapid rise of the accretion rate above the Bondi 
value, and the mean accretion rate is rather small, 
0.3Mb. The rapid and large fluctuations of the 
accretion rate appear however, and temporarily 
the accretion rate flare can reach even 0.75 of the 
Bondi rate. These changes are are caused by the 
oscillations at the outer edge of the torus, similarly 
to those found in Moscibrodzka & Proga (2008) in 
their 2-dimensional computations, however in our 
3D case the amplitude of these oscillations is much 
larger. 



In Figure [T7| we plot the tilt and twist angles 
for the model Av. The initial tilt is 0°, because 
the flow rotated initially with a positive azimuthal 
velocity. Than, the tilt increased, mostly in the in- 
nermost radii, and in the end of our simulation it 
reached 40° at 0.2 Rb. The tilt is initially differ- 
ential, and smoothly decreases with radius. The 
tilt angle for the innermost part, i.e. averaged 
from r = Tin to r = 0.05i?B, initially was periodi- 
cally changing. Despite large values of the tilt, it 
was due to the very small components of and 
Ly. After time t' > 10^, the torus started pre- 
cession, and the twist angle decreased from ~ 20° 
to ~ —20° (innermost torus) and from ~ 0° to 
~ —30° (total). We estimated that the preces- 
sion period is below ~ 6.3 x 10^, and the preces- 
sion started from f = 1.05 x 10"^. Therefore the 
trend of the tilt moment decreasing and the pre- 
cession period increasing with 7 is confirmed also 
for 7 = 1.01, if the torus is supported by the con- 
stant inflow of the rotating material from the outer 
boundary. 

3.3. Evolution of the flow for various gas 
temperatures 

We calculated the models with various posi- 
tions of the Bondi radius with respect to the 
Schwarzschild radius, which is reflected by a differ- 
ent sound speed at infinity and gas temperature. 
In the models with smaller Rb/Rs^ the evolution 
proceeds faster and the torus forms much earlier: 
at t =4x10^ for Rb/Rs = 300 and at = 2 x 10^ 
for Rb/Rs = 100. 

During the torus evolution, in these models we 
detect neither tilt nor precession. The accretion 
rate onto the black hole is much below the Bondi 



rate and it is fluctuating (Figure 18). These oscil- 
lations are caused by a strong outflow of material 
that occurs periodically because the gas wants to 
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remove some angular momentum (note that the Iq 
parameter, i.e. the circularisation radius, is the 
same in ah models). 

Nevertheless, the outflows occur symmetrically 
in these models, and in the end of the simulation 
there are no systematic asymmetries with respect 
to the equatorial plane nor the azimuth. Therefore 
for the models with large sound speed at infinity, 
the outfiow stabilizes the system with respect to 
the azimuthal perturbations and there is no pre- 
cession. The exemplary maps of the velocity field 
at the end of the simulation for the models E, F, 



G and are shown in Figure 19 



4. Discussion 

In this work, we considered the hydrodynami- 
cal models of slightly rotating, non-axisymmetric 
accretion fiows. We studied the role of adiabatic 
indices, and we verified the results from our pre- 
vious work considering the precession of the inner 
torus. 

We found that the tilt and precession occur 
not only for the gas pressure dominated fiow, de- 
scribed with the equation of state with 7 = 5/3, 
but also for the radiation dominated gas, mod- 
eled with 7 = 4/3, as wll as for smaller 7 = 1.2. 
The moment, when the torus gets tilted and starts 
precessing, correlates with 7. The period of pre- 
cession also depends on the adiabatic index, and 
is shorter for smaller 7. 

We tested also the nearly isothermal gas, with 
7 = 1.01. In this case, we could not confirm 
whether the torus precesses, in the model when 
there is no further angular momentum supply from 
the outer boundary, except for the initial condi- 
tion. For small 7, the torus is very dense and com- 
pact, and its evolution is determined by the large 
scale behaviour of the fiow. If the angular momen- 
tum is added during the simulation to the outer 
boundary of the simulation domain, such a torus 
may survive and keep the accretion rate onto black 
hole constantly smaller than the Bondi rate. Oth- 
erwise, there is a lot of material that quickly falls 
radially onto the center, and because the speed 
of sound is very small, most of this material falls 
in supersonically. The torus shrinks substantially, 
and the net accretion rate onto the black hole is 
again on the order of the Bondi rate. The accre- 
tion rate fiuctuates about this value very rapidly. 



because the shrunken torus still provides an ob- 
stacle for the radially infalling gas. The boundary 
condition that we used for most of the models in 
this work, made such a scenario of the fiow evolu- 
tion proceed very fast for small 7, and much slower 
for larger 7. 

The net outfiow of material from the equatorial 
plane is not observed in the model with 7 = 1.01, 
and at large scales the fiow is spherically symmet- 
ric (cf. Moscibrodzka & Proga 2008). However 
some tilt is visible in the innermost region, it is 
twice smaller than for 7 = 5/3 or 7 = 4/3. The 
twist angle is negligibly small, as both the and 
Ly components of the angular momentum vector 
are tiny. This angle fiuctuates periodically, show- 
ing no clear signs of precession towards any direc- 
tion. 

However, if the angular momentum is con- 
stantly added at the outer boundary, the torus 
remains large and does not shrink. In such a 
model, the tilt angle increases after t' > 10^ (units 
of dynamical times at the inner radius), and the 
torus starts precessing. Therefore, we conclude 
that the non-axisymmetric perturbations occur for 
all types of the equation of state and only the mo- 
ment when instability starts growing depends on 
7. If for small 7 the amount of infiow of the rotat- 
ing material is not sufficient, the torus may shrink 
before the effects of precession can be observed. 

The present work is a follow-up of Janiuk et al. 
(2008), where we first identified the precession of 
a torus, developing in the 3-D model of a slowly 
rotating accretion fiow. We found, that if a suf- 
ficient amount of angular momentum is added to 
the initially spherically accreting Bondi fiow, the 
instabilities rise in the innermost gas. The initial 
perturbations appear due to some small asymme- 
tries in the angular momentum distribution, which 
in the calculations can be of numerical origins, but 
to which any physical system is subject. The fiow 
then enters a transient phase of growing acous- 
tic oscillations, which manifest themselves in the 
fiuctuating sonic surface. Subsequently, the tran- 
sition to a progressive departure from the initial 
state takes place (i.e. the growing tilt of the inner 
torus and its precession about the vertical axis). 

In Janiuk et al. (2008), as well as in the present 
work, we checked that the instability develops in 
a highly supersonic fiow, whereas for small Mach 
numbers it is suppressed. We can attribute this 
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behaviour with the Papaloizou & Pringle (1985) 
type of instabihty, which for low order modes is 
driven by the Kelvin-Helmholtz mechanism. How- 
ever here, the picture is more comphcated due to 
the supersonic nature of the flow, presence of out- 
flows and shocks, as well as compressibility. The 
present work was devoted to the effect of the var- 
ious adiabatic index 7. 

We found that the inner torus precesses also for 
values other than 7 = 5/3, used in our first study: 
7 = 4/3, 1.2, and 1.01. However, the time for the 
precession to set increases with decreasing 7. For 
the nearly isothermal model, 7 = 1.01, the time of 
the instability growth was too long to set the pre- 
cession, unless the outer boundary condition was 
changed and more angular momentum was sup- 
plied during the time evolution. Our results are 
in agreement with the general behaviour of the 
unstable accreting flows for various equations of 
state. Blondin et al. (2003) studied the stability of 
standing accretion shocks, and also found that the 
softer the equation of state (smaller 7), the slower 
is the growth of SAS instability. For the smallest 
7 = 1.25 used by these authors, they found that 
the shock remained marginally stable for a long 
period of time, however also this model eventually 
became as unstable as the others. The physical 
mechanism that is lying behind this SASI insta- 
bility is the vortical-acoustic feedback between the 
pressure waves and vorticities induced by the as- 
pherical shock. Such feedback was also studied by 
Foglizzo (2002) in the context of accreting black 
holes. In our simulations, shocks arise due to the 
outflow of gas in the equatorial plane, and no po- 
lar outflows are produced without the magnetic 
field (see e.g. Moscibrodzka & Proga 2009), but 
the mechanism appears to be similar. 

The mechanism of precession that we present in 
this work appears to be an additional one, apart 
from those discussed in the litearture. Caproni 
et al. (2006) discuss several different sources of 
the accretion disk warping and precession. It may 
occur due to the tidal forces, irradiation, magnet- 
ically driven instability, or the Bardeen-Petterson 
effect. Signatures of the disk precession are ob- 
served in both AGN and stellar mass black hole 
or neutron star binaries. The observations involve 
the spectral and temporal variability as well as 
the distortions of jet morphology. The period of 
precession can be quite short, and for the sources 



with black hole mass on the order of 10^ - lO^M© 
the observed times are on the order of 1-10 years. 
Our results are in good agreement with these ob- 
servations of AGN. 

To study the precession mechanism in more de- 
tail we need to include the magnetic fields in our 
modeling. If the Kelvin-Helmholtz mechanism lies 
behind the torus precession instability, the results 
should be affected by the presence of magnetic 
field component only in the streaming direction. 
This is planned to be the subject of our future 
work. In addition, it will be useful to compare 
the results of the present ZEUS-MP code with an- 
other MHD code available for testing, to help clar- 
ify whether the main instability in our studies is 
of a physical or numerical nature. 
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Fig. 2. — Time evolution of the angular momen- 
tum flux through the inner boundary, in units of 
the critical angular momentum times the Bondi 
accretion rate. The models are: A (red line), B 
(green line), C (blackline), D (cyan line). 
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Fig. 3. — Contour density maps and velocity field 
in the inner region for the model A (7 = 1.01), in 
the side-view (along the x-axis; left) and in the top 
view (along the z-axis; right). The time snapshots 
are taken at = 3 x 10^,6.3 x 10^ and 1.1 x 10"^, 
from top to bottom. The units on the x, y and z 
axes are in Rb- 



Fig. 4. — Contour plots of the total Mach num- 
ber, Mtot (left) and the radial Mach number, 
(right) for the model A (7 = 1.01). The time 
snapshots are taken at t = 3 x 10^, 6.3 x 10^ and 
1.1 X 10^, from top to bottom. The units on the 
X and z axes are in Rb, (the bottom-right panel 
is zoomed-in). The color contours show the con- 
stant Mach number of 0.5, 1.0 and 1.5 (red, green 
and blue, respectively), while the grey colors show 
the constant density contours, between 10^ and 

lO^oo. 
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Fig. 5. — Contour density maps and velocity field 
in the inner region for the model B (7 = 1.2), in 
the side-view (along the x-axis; left) and in the top 
view (along the z-axis; right). The time snapshots 
are taken at = 3 x 10^,6.3 x 10^ and 1.1 x 10"^, 
from top to bottom. The units on the x, y and z 
axes are in Rb- 



Fig. 6. — Contour plots of the total Mach number, 
M (left) and the radial Mach number, Mr (right) 
for the model B (7 = 1.2). The time snapshots 
are taken at = 3 x 10^,6.3 x 10^ and 1.1 x 10"^, 
from top to bottom. The units on the x and z axes 
are in Rb. The color contours show the constant 
Mach number of 0.5, 1.0 and 1.5 (red, green and 
blue, respectively), while the grey colors show the 
constant denity contours, between 10^ and lO^Poo- 
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Fig. 7. — Contour density maps and velocity field 
in the inner region for the model C (7 = 4/3), in 
the side-view (along the x-axis; left) and in the top 
view (along the z-axis; right). The time snapshots 
are taken at = 3 x 10^,6.3 x 10^ and 1.1 x 10"^, 
from top to bottom. The units on the x, y and z 
axes are in Rb- 



Fig. 8. — Contour plots of the total Mach number, 
M (left) and the radial Mach number, Mr (right) 
for the model C (7 = 4/3). The time snapshots 
are taken at = 3 x 10^,6.3 x 10^ and 1.1 x 10"^, 
from top to bottom. The units on the x and z axes 
are in Rb. The color contours show the constant 
Mach number of 0.5, 1.0 and 1.5 (red, green and 
blue, respectively), while the grey colors show the 
constant denity contours, between 10^ and lO^Poo- 
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Fig. 9. — Contour density maps and velocity field 
in the inner region for the model D (7 = 5/3), in 
the side-view (along the x-axis; left) and in the top 
view (along the z-axis; right). The time snapshots 
are taken at = 3 x 10^,6.3 x 10^ and 1.1 x 10"^, 
from top to bottom. The units on the x, y and z 
axes are in Rb- 



Fig. 10. — Contour plots of the total Mach num- 
ber, M (left) and the radial Mach number, 
(right) for the model D (7 = 5/3). The time 
snapshots are taken at t = 3 x 10^, 6.3 x 10^ and 
1.1 X 10^, from top to bottom. The units on the 
X and z axes are in Rb- The color contours show 
the constant Mach number of 0.5, 1.0 and 1.5 (red, 
green and blue, respectively), while the grey col- 
ors show the constant denity contours, between 
10^ and lO^oo. 
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Fig. 11. — Sequences of the poloidal velocity 
fields, taken for the azimuth (/> = 0, for the mod- 
els with 7 = 1.01,1.2,4/3 and 5/3 (from left 
to right) and at three different time snapshots: 
t = 3 X 10^6.3 X 10^ and 1.1 x 10^ (from top 
to bottom). 



Fig. 12. — The tilt (left panel) and averaged twist 
(right panel) angles, calculated for model with 7 = 
5/3. The tilt is shown as a function of radius, 
for several time-snapshots, given in the top of the 
panel. The initial tilt is of 180°, meaning that the 
angular momentum vector is not tilted to the z- 
axis and the flow rotates with a negative azimuthal 
velocity. The twist angle is shown as a function of 
time, and the two curves represent the averaging 
over radius from r = rin to r = 0.05i?B (solid 
line) and to r = 1.2Rb (dashed line). The twist 
angle decreasing from ~ 60° to ~ —80° (innermost 
torus) and from ^ 40° to ^ -40° (total) reflects 
the clockwise, initially differential precession. 
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Fig. 13. — The tilt (left panel) and averaged twist 
(right panel) angles, calculated for the model with 
7 = 4/3. The tilt is shown as a function of radius, 
for several time-snapshots, given in the top of the 
panel. The initial tilt is of 180° means that the 
angular momentum vector is not tilted to the z- 
axis and the flow rotates with a negative azimuthal 
velocity. The twist angle is shown as a function of 
time, and the two curves represent the averaging 
over radius from r = rin to r = 0.05i?B (solid line) 
and to r = 1.2i?B (dashed line). The twist angle 
increasing from ~ —120° to ~ —40° (innermost 
torus) and from ~ —70° to ~ 0° (total) reflects the 
counter-clockwise, initially differential precession. 
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Fig. 14. — The tilt (left panel) and averaged twist 
(right panel) angles, calculated for the model with 
7 = 1.2. The tilt is shown as a function of radius, 
for several time-snapshots, given in the top of the 
panel. The initial tilt is of 180° means that the an- 
gular momentum vector is not tilted to the z-axis 
and the flow rotates with a negative azimuthal ve- 
locity. The twist angle is shown as a function of 
time, and the two curves represent the averaging 
over radius from r = rin to r = 0.05i?B (solid line) 
and to r = 1.2Rb (dashed line). The twist angle 
increasing from ~ —150° to ~ —10° (innermost 
torus) and from ~ —80° to ^ 10° (total) reflects 
the counter-clockwise, initially differential preces- 
sion. 
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Fig. 15. — The tilt (left panel) and averaged twist 
(right panel) angles, calculated from the time evo- 
lution of the model A (7 = 1.01). The tilt is shown 
as a function of radius, for several time-snapshots, 
given in the top of the panel. The initial tilt is 
of 180° means that the angular momentum vector 
is not tilted and the flow rotates with a negative 
azimuthal velocity. The twist angle is shown as a 
function of time, and the two curves represent the 
averaging over radius from r = rin to r = 0.05i?B 
(solid line) and to r = 1.2Rb (dashed line). The 
angular momentum vector is always almost per- 
pendicular to the equatorial plane, and despite the 
large scatter in the two angles, the apparent small 
tilt is due to the very small components of and 

Ly. 
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Fig. 16. — Time evolution of the mass accretion 
rate through the inner boundary (Min), in units 
of the Bondi rate (Mb). The model is Av^ i.e. 
7 = 1.01 but with the changed outer boundary 
condition. 




■■ [R.] time [V(r,.)] 



Fig. 17. — The tilt (left panel) and averaged twist 
(right panel) angles, calculated from the time evo- 
lution of the model Av (7 = 1.01, but with the 
changed outer boundary condition). The tilt is 
shown as a function of radius, for several time- 
snapshots, given in the top of the panel. The 
initial tilt of 0° means that the angular momen- 
tum vector is not tilted and the flow rotates with 
a positive azimuthal velocity. The twist angle is 
shown as a function of time, and the two curves 
represent the averaging over radius from r = rin 
to r = 0.05i?B (solid line) and to r = 1.2Rb 
(dashed line). The twist angle decreasing from 
- 20° to - -20° (innermost torus) and from - 0° 
to ~ —30° (total) reflects the clockwise, initially 
differential precession, starting from time t' > 10^. 
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Fig. 18. — Time evolution of the mass accretion 
rate through the inner boundary (Min), in units 
of the Bondi rate (Mb). The models are: E 
(red line), F (green line), G (black line) and H 
(cyan line). The models have the adiabatic in- 
dex of 7 = 5/3 and 7 = 4/3, as labelled on the 
plot, and temperature determined by the ratio of 
R' = Rb= Rs is equal to 100 or 300. 
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Fig. 19. — The poloidal velocity field for the (j) = 
slice, at the end of the simulation. The models are: 
E (top left), F (top right), G (bottom left) and H 
(bottom right). 
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